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Abstract: Fracture processes within a material begin at the nanometer length 
scale at which the formation, propagation, and interaction of fundamental damage 
mechanisms occur. Physics-based modeling of these atomic processes quickly 
becomes computationally intractable as the system size increases. Thus, a 
multiscale modeling method, based on the aggregation of fundamental damage 
processes occurring at the nanoscale within a cohesive zone model, is under 
development and will enable computationally feasible and physically meaningful 
microscale fracture simulation in polycrystalline metals. This method employs 
atomistic simulation to provide an optimization loop with an initial prediction of a 
cohesive zone model (CZM). This initial CZM is then applied at the crack front 
region within a finite element model. The optimization procedure iterates upon 
the CZM until the finite element model acceptably reproduces the near-crack- 
front displacement fields obtained from experimental observation. With this 
approach, a comparison can be made between the original CZM predicted by 
atomistic simulation and the converged CZM that is based on experimental 
observation. Comparison of the two CZMs gives insight into how atomistic 
simulation ‘scales.’ 

1 Introduction 

State-of-practice methods for predicting material fracture are based on empirical 
models that homogenize damaging mechanisms and recast them in the form of a 
single fracture or damage parameter, e.g, Crack Tip Opening Angle (CTOA). In 
general, a damage parameter is asserted to be the driving force for crack 
propagation, and subsequently, mechanical testing is performed to provide the 
corresponding critical value for a given material and testing condition. State-of- 
practice methods for modeling fracture have greatly improved the ability to 
predict structural failure; however, these methods are inductive and only a small 
subset of all possible configurations (both material and geometric) can be tested. 
This results in high uncertainty and safety factors, and necessitates decreased 
inspection intervals and inefficient design. The state-of-practice methods are 
fundamentally limited in that a phenomenological governing driving force, such 
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as CTOA, is initially asserted, and testing is carried out to find the corresponding 
value of critical material resistance, i.e. critical CTOA (CTOA c ), that was initially 
asserted to be the mechanism of consequence. Breaking out of this circular- 
reasoning loop and developing a more deductive approach enabled by modeling 
fundamental mechanisms will reduce uncertainty and result in more efficient and 
reliable structures. 

Useful prognosis of structural behavior is probabilistic due to the inherent 
variability in materials, loading, manufacturing, and environment at all physical 
scales. In a statistical sense, lower uncertainty in prognoses implies that this 
inherent variability is understood and incorporated within the analysis. It is the 
contention of the authors that a statistical prognosis capability requires 
quantitative simulation of microstrutural crack growth utilizing physics-based 
models of intra-, inter-, and transgranular fracture modes for polyciystalline 
metals. However, it is computationally intractable and arguably unnecessary to 
represent an entire component with atomistic resolution, especially statistically. 
Therefore, an approach that leads to mechanistically-based continuum fracture 
parameters developed via multiscale modeling and various forms of 
homogenization is desired. 

In this paper, a methodology for development of a physically-grounded analysis 
for modeling the initial stages of crack growth in an aluminum alloy is described. 
Section 2 of this paper discusses a method for improvement of CZMs that is based 
on the aggregation of experimentally-derived and computationally-simulated 
response of single- and bi-crystals. In Section 3, a comparison between state-of- 
practice methods and the methodology presented in Section 2 is given. Section 3 
serves as an exercise to illustrate the expected advances to be made through the 
incorporation of nanoscale mechanics (discontinuum) into microscale simulation 
(continuum). 

2 Molecular dynamics, experimentation, and inverse modeling 

The purpose of molecular dynamics (MD) simulation in this study is to reveal and 
analyze nanoscale processes occurring near a crack front so that they can be 
incorporated into crack growth simulation through CZMs. As shown in past 
work, MD simulation can be used to determine a CZM in a statistical sense 
[8,11,12]. From MD simulation, one can determine the influence that various 
plastic processes have on crack growth, such as vacancy and void formation, and 
dislocation nucleation. Crack growth simulation via a finite element model 
(FEM) can then be improved by incorporating a CZM that takes into account 
these distinct processes near the crack front. Currently, however, MD simulation 
is applied to idealized crystal structures composed of materials for which 
interatomic potentials have been developed and, as such, the analyses tend to be 
of only qualitative use. Also, current compute power limits the simulations to 
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very short time scales (fractions of a microsecond) and small spatial domains 
(fractions of a micron) [10]. 

Fig. 1 illustrates the extraction of a MD-based CZM representing decohesion 
along an idealized E99 grain boundary in an aluminum bi-ciystal [11]. While 
these MD-predicted CZMs are qualitatively similar to experimentally derived 
CZMs, they are not quantitatively similar. This is a result of the idealized nature 
of MD simulation. From Fig. 1, it is seen that MD simulation predicts a CZM 
with an ultimate tensile traction, ct c , of nearly 6 GPa and a complete separation 
distance, X c , of approximately 2.4 nm [11] for pure aluminum, whereas in 
practice, values of around 500 MPa and 20 pm are expected for a c and X c , 
respectively, for aluminum alloys [2,4], The ability to relate a CZM as predicted 
by MD simulation to a CZM observed through physical experimentation is 
needed to develop a validated model of the process zone of a crack; such a 
capability would provide insight into how molecular dynamics simulation 
‘scales.’ 



Fig. 1: CZM extraction from MD simulation [11] 

To best match the MD simulation spatially, small-scale experiments and highly 
resolved measurement techniques must be employed. Fig. 2 is an illustration of 
such an experiment where a fatigue-cracked, single-crystal of 6019 aluminum 
alloy is being loaded uniaxially in the y-direction, and strains in the y-direction 
are measured using high-resolution (approximately 10 nm displacement 
resolution) digital image correlation (DIC). In the experimental approach 
advocated here, single-crystals of pure aluminum, rather than an aluminum alloy, 
are fabricated and fatigue-cracks are initiated. Upon loading, displacements are 
measured using DIC in the area surrounding the point at which the crack front 
intersect the specimen surface. This near-crack displacement field enables 
observation of the localization of strain at the crack front. MD modeling 
techniques can then be employed to determine the role that each mechanism plays 
in the process zone and subsequently its contribution to the CZM. 

While MD simulation is capable of providing qualitatively accurate CZMs that 
account for the contribution of distinct mechanisms, and experimentation can 
provide high-resolution displacement fields, there is no reconciliation of the 
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information gained. This reconciliation is the key step to developing a validated 
model for the process zone. Unfortunately, it is currently computationally 
intractable to extend the temporal and spatial domains of MD simulation beyond 
the environs discussed previously, and experimentation certainly cannot measure 
the movement of each atom, especially on a picosecond scale. Therefore, the 
methodology developed here aggregates the information gained through MD 
simulation for application at scales, both temporal and spatial, which lend 
themselves to experimental validation. 



Fig. 2: Digital image correlation produced contours ofy-direction strain 

At this point, an inverse problem has manifested itself; the displacement field has 
been measured, but one does not know which material model produces that 
displacement field. Past research regarding the inverse calculation of CZMs to 
simulate these forms of crack growth has focused, to a large extent, on numerical 
methods for solving the resulting nonlinear systems [3]. This inevitably has led to 
a choice of CZM parameters that is numerically advantageous rather than being 
physically meaningful, since the mechanisms that control those parameters are not 
well-known and have not been quantified. Concrete and ceramics have been the 
focus of much past research on calculating CZM parameters using inverse 
methods. Previously developed methods for these materials are macroscopic in 
nature, in that the CZM functional form is fixed (usually piecewise linear) and its 
parameters are calculated using far-field displacement measurements 
[1,2, 4, 6, 7, 9]. Using these macroscopic approaches for inverse calculation of 
CZMs reproduces well the far-field displacement and total energy dissipation of a 
system. However, the major goal of the present work is to determine the impact 
that distinct processes, such as those observed by Horton [5] and others, have on 
the CZM. 

A technique for the inverse calculation of a CZM, Fig. 3, is now formulated using 
the experimentally measured displacement field near the crack process zone and 
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an initial guess for the CZM, based on the MD simulation. These initially 
predicted CZMs are placed within a finite element model that attempts to replicate 
the configuration, continuum material response, boundary conditions, and loading 
seen in the experiment. A field of near-crack, surface displacements is measured 
using the experimental setup and used as the reference solution, u e . The loop, 
denoted by the circle in Fig. 3, continues until the corresponding finite element 
displacement field, u s , is within a defined tolerance of the experimental field, u e . 
The inverse problem-solving technique converges to a CZM that produces the 
best replication of experimental observation. 


The optimization technique employed here for development of the CZM is based 
on a genetic algorithm. The genetic algorithm technique is employed because of 
its inherent capability to search for global minima and its ease of parallelization. 
Initially, a set of CZMs is generated that is closely related to the MD-predicted 
CZM. An additional set is also generated that varies the CZM parameters by an 
order of magnitude, which enables a global search of optimal CZM parameters. 
For each CZM in the population, a finite element model is run; using a computer 
cluster, each FEM is sent to a different processor so that runtime is greatly 
reduced. Once the comparison is made between each CZM in the population and 
the experimentally observed field, if convergence has not been achieved, the 
genetic algorithm is called to update the population of CZM parameters. It is 
understood that the outcome of this process is not unique and that the converged 
solution can, for ill-posed problems, simply be a function of the initial guess. 
Therefore, subsequent studies will be carried out that investigate the convexity of 
the search space and the sensitivity to initial predictions of the CZM. 
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Fig. 3: Inverse approach to finding CZM parameters 
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3 An illustrative example 


In this section, a comparison between state-of-practice methods and the 
methodology under development is provided. A bicrystal model undergoing 
intergranular crack growth is used as the basis for comparison, Fig. 4, and was 
chosen to illustrate the advantage of using a MD-informed CZM during crack 
growth simulation. Upon loading uniaxially in the y-direction, the two crack tips 
will begin developing process zones. Let the orientations of Crystal #1 and #2 be 
such that the interface is symmetric, as in [11,12], The goal of this comparison is 
to illustrate the conceptual differences amongst modeling approaches for crack 
growth simulation rather than to detail quantitative comparison or validation. 


In the analogous FEM simulation, the two crystals are assumed to be elasto- 
plastic with power-law hardening, i.e. a Ramberg-Osgood relation, so that stress 
and strain are given by the Hutchinson, Rice, and Rosengren (HRR) fields, Eqs. 1. 
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Eqs. 1 


In Eqs. 1: 

Ojj = stress tensor 
= strain tensor 

< % , I„ = dimensionless constants 

o n = yield stress 

n = strain hardening exponent 

0 = angle about crack tip, Fig. 4 

dy (0,n) ,£-(0,11) = dimensionless functions of 9 and n 

The state-of-practice methods are well-represented by three different approaches 
to this problem: 

Damage mechanics - where a critical region or distance is defined from the crack 
front and a computed damage parameter is compared to a measured critical value 
to determine whether or not the crack will grow and; 

Crack opening displacement (COD) - instead of querying a damage parameter at 
a distance ahead of the crack, opening is queried at a distance behind the crack tip 
(conceptually not much different from damage mechanics) and; 
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Cohesive zone modeling (CZM) - a planar material model is placed at the crack 
front to account for plastic processes localized there, the defining parameters are 
commonly chosen because they are numerically advantageous or are the result of 
experimentation or an inverse technique [3], 



Fig. 4: Idealized bicrystal undergoing intergranular crack growth 

Substitution of the bounding values of n into Eqs. 1 reveals some fundamental 
limitations of state-of-practice methods. First, if n = 1, it is found that both stress 
and strain are proportional to l/(r) 1/2 (i.e. the linear-elastic case). Conversely, as n 
-> oo, then strain is proportional to 1/r and stress is effectively independent of r 
{i.e. the perfectly-plastic case). At either bound, strain fields are singular, which 
is directly a result of defining an infinitely sharp crack in the derivation of the 
HRR fields. This singularity is physically invalid in an elasto-plastic material due 
to crack tip blunting, which means that the HRR stress and strain fields are not 
valid at the crack front, r = 0. This leads to the motivation for both damage 
mechanics and COD. Since HRR fields are invalid at the point where they are 
most important, one can move some distance away and query fields there (or 
average over some region). 

However, these methods inevitably lead to the need for nonphysical parameters, 
such as a critical distance, and therefore need to be calibrated for each geometrical 
and loading configuration considered. In addition to the use of nonphysical 
parameters, a problem exists with mesh convergence. Since strain is singular at 
the crack tip, and the form of the singularity depends on Q , n, grain orientation, 
etc., a converged solution at the crack front is not possible with mesh refinement. 
So, another calibration must occur where the finite element size is a fixed ratio of 
the critical distance, providing a basis for comparison among various finite 
element models. However, for complicated geometry, (as is the case for 
microstructural geometry) variations in mesh refinement are needed to obtain 
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reasonable element shapes and, in some cases, even to make meshing possible. 
So, although this is a logical approach when comparing idealized models within a 
well-defined study, it is not a general solution, especially for complicated 
microstructural geometry. 

Crack growth simulation based on CZMs does not produce a singularity at the 
crack tip as occurs in damage mechanics or COD-based methods. As discussed in 
Sec. 2, several recent publications detail methods for back-calculation of CZM 
parameters from far-field displacements. However, in the case shown in Fig. 4, 
measurement of far-field displacements does not help to distinguish among the 
various physical processes at each tip and necessarily leads to a CZM for both 
process zones that does not account for the dominant mechanisms local to the 
particular crack front. The mechanisms for plasticity within the process zones 
are, in fact, very different at each tip. Yamakov, et al., [11,12] showed, for an 
embedded crack growing along a S99 grain boundary, that one tip will grow in a 
relatively brittle manner, through void formation, while the other will grow in a 
more ductile manner, with twin emission blunting the crack tip. This hindrance 
can be overcome by employing the experimental and measurement methods 
presented in Sec. 2. By directly measuring the displacement field about each 
process zone, the back-computed CZMs can be separated and more accurately 
depict the physical processes occurring at each tip. Not only does the initial MD 
simulation of such a problem make the understanding of various dominant 
mechanisms more obvious, it also allows for characterization of process zones 
through CZMs by the mechanisms that control the local deformation. 

4 Conclusions 

Given the need for a more physics-based approach to computationally modeling 
crack growth, there is strong motivation to begin incorporating MD simulation as 
part of a multiscale analysis. However, current MD simulations are based on 
unrealistic idealizations; simulations occur over such small length and time scales 
that they cannot be directly compared with experiment. On the other hand, MD 
simulation provides useful qualitative insight into the physical processes 
occurring at a crack. This insight can be used to better predict and understand 
CZMs modeling the process zone. 

Past research has also motivated refinements in experimental observations and 
measurement. Combining this refinement in measurement techniques with 
corresponding numerical methods can provide near-crack displacement fields and 
more accurate models of the process zone through CZMs. In addition, with an 
understanding of the distinct mechanisms occurring within the process zone from 
MD simulation, the possibility to validate MD simulation exists. Development of 
CZMs in this way will improve the existing capability to model the early stages of 
crack growth. 
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